function [CHECKvariable] = extractModels(subfolder,SpreadsheetName,TabName_Step2,WorksheetName2Save)

CHECKvariable=0;
mkdir(subfolder);
addpath(subfolder)


warning('off','MATLAB:MKDIR:DirectoryExists')



%% Step 1: import the data from the spreadsheet, simultaneously normalize and calculate new parameters (such as Mg#).
addpath(genpath(pwd));

%[type, TabNames] = xlsfinfo(SpreadsheetName);
Labels2Plot = [];
TabNames={...
    'MASTER'...
    };


PetrogenMajorElement_Strings = {'SiO2','TiO2','Al2O3','Cr2O3','FeO','MgO','CaO','K2O','Na2O'};
% PetrogenTraceElement_Strings
% Petrogen_Co
% Petrogen_Co_Name
% Petrogen_BulkCompositionsNormalized
% Petrogen_Dmatrix

targetStrings_Majors = {'T(C)', 'P(kbars)','SiO2','TiO2','Al2O3', 'Cr2O3','FeO','MnO','MgO','CaO','Na2O','K2O','P2O5','NiO','H2O'};

ElementRatios_Strings= {'Th' 'U'; 'Lu' 'Hf'; 'Sm' 'Nd'; 'Lu' 'Sc'; 'Hf' 'Nd'; 'La' 'Sm'; 'Rb' 'Sr'; 'Eu' 'Sm'; 'Eu' 'Gd'; 'Nd' 'Zr'; 'Sm' 'Yb'; 'Zr' 'Y';...
    'Nb' 'Y'; 'Gd' 'Yb';  'Zr' 'Yb'; 'Zr' 'Hf'; 'Nb' 'Ta'; 'Ba' 'La'; 'Ba' 'Pb'; 'U' 'Pb'; 'Eu' 'Eu'; 'K' 'U'};

% MajorElements_Strings = {'T(C)' 'P(kbars)' ...
%     'SiO2'	'TiO2'	  'Al2O3' 'Cr2O3'	  'FeO'	  'MnO'	  'MgO'	  'CaO' 'Na2O'   'K2O' 'P2O5' 'NiO' 'H2O' 'total'...
%     'Mg#' 'NaK#' 'CaO/Al2O3 wt%' 'CaO/Al2O3 moles' 'Na2O/FeO' 'K2O/TiO2' '1-Mg#' 'Qtz' 'Plag' 'Oliv' 'Cpx' 'Ox' 'Or' 'Ap'};



% Loads in constraint option names
[xlsNumbers, xlsText,xlsRAW] = xlsread(SpreadsheetName,'HeadingOptions');
[ReferenceRow,ReferenceColumn]= find(strcmp(xlsText,'FirstHeading')==1);
[LastReferenceRow,LastReferenceColumn]= find(strcmp(xlsText,'LastHeading')==1);
ConstraintOptions = xlsText(ReferenceRow+1, ReferenceColumn:LastReferenceColumn);
% ConstraintOptions = {'dFdPwaterchange','rFrac_crit','initialwatercontent','Tmp','U','BulkComp',...
%     'FirstMeltPb','MeltExtractedPb','LastMeltPb','Gar2SpTransition','Sp2PlagTransition',...
%     'aluminumOutPb','aluminumOutFrac','cpxOutPb','cpxOutFrac','Tdiff','quick crustal thickness',...
%     'ct2_ininc(1)','percentGar','percentSpinel','percentPlag','Full Pool Fmean','TYPE',...
%     'percentPool','distance_ininc','Variable Pool Fmean','Vpool %gar','Vpool  %sp','Vpool %plag','gar km','sp km','plag km','totalkm',...
%     'Pb','T','Ts','Ln_cumulative','Lnexp_cumulative','MeltType'}



%Loads in target headings for this extraction
[ReferenceRow,ReferenceColumn]= find(strcmp(xlsText,'FirstExtractHeading')==1);
[LastReferenceRow,LastReferenceColumn]= find(strcmp(xlsText,'LastExtractHeading')==1);
Info_Order = xlsText(ReferenceRow+1, ReferenceColumn:LastReferenceColumn);
% THIS IS WHAT YOU WILL EXTRACT:
% Info_Order = {'dFdPwaterchange','rFrac_crit','Water','Tmp','U','Bulk Comp','FirstMeltPb','MeltExtractedPb','LastMeltPb','Gar2SpTransition','Sp2PlagTransition',...
%     'percentGar','percentSpin','percentPlag','ct2_ininc(1)',...
%     'Full Pool Fmean','Full Pool Fmean','Tdiff','aluminumOutPb','aluminumOutPb','aluminumOutFrac','cpxOutPb','cpxOutPb','cpxOutFrac'};


% 'OnAxis'...
% 'FullPool'...
% 'GarPool'...
% 'FirstMelt'...


%% Step 0: spectify your target major and trace element orders

%%



[xlsNumbers, xlsText,xlsRAW] = xlsread(SpreadsheetName,TabName_Step2);
[ReferenceRow,ReferenceColumn]= find(strcmp(xlsText,'FirstTargetName')==1);
[LastReferenceRow,LastReferenceColumn]= find(strcmp(xlsText,'LastTargetName')==1);
ModelNames = xlsText(ReferenceRow+1:LastReferenceRow-1, ReferenceColumn);
numberReferences = size(ModelNames,1);

[ElementRow,ElementColumn]= find(strcmp(xlsRAW,'FirstTarget')==1);
[LastElementRow,LastElementColumn]= find(strcmp(xlsRAW,'LastTarget')==1);
TargetList = xlsRAW(ElementRow+1, ElementColumn:LastElementColumn);
TargetListData = xlsRAW(ElementRow+2:LastReferenceRow-1, ElementColumn:LastElementColumn);
TargetListData=cell2mat(TargetListData);

[Step1Row, Step1Column]= find(strcmp(xlsText,'Step1MatFile')==1);
NeededStep1Files  = (xlsText(ReferenceRow+1:LastReferenceRow-1, Step1Column));
NeededStep1Files_unique = unique(NeededStep1Files);


[a,infoIndicies]=ismember(Info_Order, ConstraintOptions);

   TargetList_table = regexprep(TargetList,' ','');
   OutputTable = table(ModelNames,NeededStep1Files,TargetListData);
   OutputTable = splitvars(OutputTable,'TargetListData', 'NewVariableNames',TargetList_table);
   
   disp(OutputTable)
   checkInput = input(sprintf('\n   <strong>If these are wrong, enter ''N'', otherwise any key to continue:</strong>',TabName_Step2),'s');
   
   if strcmp(checkInput,'N')==1
       CHECKvariable = 1; 
       return
   end
   
Names2Saves ={};

All_INFO = [];

for numSteps = 1:size(NeededStep1Files_unique,1)
   
    PETROGEN_output2load = NeededStep1Files_unique{numSteps};
    uniqueStepRefs = find(ismember(NeededStep1Files, PETROGEN_output2load))';
    
    
    load(PETROGEN_output2load)
  
    %Get indicies for element ratios
    EuCol = find(strcmp(PetrogenTraceElement_Strings,'Eu'));
    SmCol = find(strcmp(PetrogenTraceElement_Strings,'Sm'));
    GdCol = find(strcmp(PetrogenTraceElement_Strings,'Gd'));
    for gg = 1:size(ElementRatios_Strings,1)
        for ggg = 1:size(ElementRatios_Strings,2)
            ElementRatiosIndex(gg,ggg) =  find(strcmp(PetrogenTraceElement_Strings,ElementRatios_Strings(gg,ggg)));
        end
    end
    
    
    for kk =uniqueStepRefs
        
        Data2Find = TargetListData(kk,:);
        ReferenceName = ModelNames{kk};
        
        for TabMelts=1:size(TabNames,2)
            
            
            
            InfoElement_ALL = eval(sprintf('%s_INFO',TabNames{TabMelts}));
            InfoNames = ConstraintOptions;
            [a,b]=ismember(TargetList,InfoNames);
            data2plotAll_TargetInfo= InfoElement_ALL(:,b);
            
          
            compIndicies=find(any(bsxfun(@minus,data2plotAll_TargetInfo,Data2Find),2)==0)';
     
            clear data2reorder ReorderedData
            data2reorder = eval(sprintf('%s_MAJORS',TabNames{TabMelts}));
            data2reorder = data2reorder(compIndicies,:);
            [A,ElementIndicies4Target] = ismember(targetStrings_Majors, PetrogenMajorElement_Strings);
            [reformattedmajors,tormeyhere,PetrogenMajorElement_Strings_Extracted]=sumMgNumNorm_PT(data2reorder,ElementIndicies4Target,[98.5 101.5]);
       
            
         
            
            
%             [A,ElementIndicies4Target] = ismember(targetStrings_Majors, PetrogenMajorElement_Strings);
%             nodata = find(ElementIndicies4Target==0);
%             ElementIndicies4Target(ElementIndicies4Target == 0) = max(ElementIndicies4Target);
%             ReorderedData = data2reorder(:,ElementIndicies4Target);
%             %pads rows of NaNs for elements with missing data
%             ReorderedData(:,nodata)=[NaN];
%             MajorElements=ReorderedData(compIndicies,:); 
%             reformattedmajors = sumMgNumNorm(MajorElements);
%             
            
           
            
            
            
            
            
            InfoElement=InfoElement_ALL(compIndicies,infoIndicies);
            TraceElements_ALL = eval(sprintf('%s_TRACE',TabNames{TabMelts}));
            TraceElements = TraceElements_ALL(compIndicies,:);
            
            
            Residue_Trace_ALL = eval(sprintf('%s_Residue_TRACE',TabNames{TabMelts}));
            Residue_Trace = Residue_Trace_ALL(compIndicies,:);
            
            SolidResidue_Trace_ALL = eval(sprintf('%s_SolidResidue_TRACE',TabNames{TabMelts}));
            SolidResidue_Trace = SolidResidue_Trace_ALL(compIndicies,:);
            
            CPX_ALL = eval(sprintf('%s_CPX',TabNames{TabMelts}));
            CPX = TraceElements_ALL(compIndicies,:);
            
            Ds_ALL = eval(sprintf('%s_Ds',TabNames{TabMelts}));
            Ds = Ds_ALL(compIndicies,:);
            
            DISEQ_ALL = eval(sprintf('%s_DISEQ',TabNames{TabMelts}));
            DISEQ = DISEQ_ALL(compIndicies,:);
            
            
            %labelhere =  sprintf('%s_%s',ReferenceName,TabName2Save);
            labelhere =  sprintf('%s',ReferenceName);
            
          
           
            clear DataRatios Residue_DataRatios SolidResidue_DataRatios CPX_DataRatios
            for gg = 1:size(ElementRatios_Strings,1)
                DataRatios(:,gg) = TraceElements(:,ElementRatiosIndex(gg,1)) ./TraceElements(:,ElementRatiosIndex(gg,2));
                DataRatios(:,21) = 2.*(TraceElements(:,EuCol)./0.154)./((TraceElements(:,SmCol)./0.406)+(TraceElements(:,GdCol)./0.544));
                
                Residue_DataRatios(:,gg) = Residue_Trace(:,ElementRatiosIndex(gg,1)) ./Residue_Trace(:,ElementRatiosIndex(gg,2));
                Residue_DataRatios(:,21) = 2.*(Residue_Trace(:,EuCol)./0.154)./((Residue_Trace(:,SmCol)./0.406)+(Residue_Trace(:,GdCol)./0.544));
                
                SolidResidue_DataRatios(:,gg) = SolidResidue_Trace(:,ElementRatiosIndex(gg,1)) ./SolidResidue_Trace(:,ElementRatiosIndex(gg,2));
                SolidResidue_DataRatios(:,21) = 2.*(SolidResidue_Trace(:,EuCol)./0.154)./((SolidResidue_Trace(:,SmCol)./0.406)+(SolidResidue_Trace(:,GdCol)./0.544));
                
                
                CPX_DataRatios(:,gg) = CPX(:,ElementRatiosIndex(gg,1)) ./CPX(:,ElementRatiosIndex(gg,2));
                CPX_DataRatios(:,21) = 2.*(CPX(:,EuCol)./0.154)./((CPX(:,SmCol)./0.406)+(CPX(:,GdCol)./0.544));
                
                TESource_DataRatios(:,gg) = Petrogen_Co(ElementRatiosIndex(gg,1)) ./Petrogen_Co(ElementRatiosIndex(gg,2));
                TESource_DataRatios(:,21) = 2.*(Petrogen_Co(EuCol)./0.154)./((Petrogen_Co(SmCol)./0.406)+(Petrogen_Co(GdCol)./0.544));
                
            end
            
            
           
            % Must call a new script to properly save the melt compositions
            savevariables
            
            
            
            % compIndicies=find(any(bsxfun(@minus,data2plotAll_TargetInfo,Data2Find),2)==0)';
            % clear data2reorder ReorderedData
            % data2reorder = eval(sprintf('%s_SM_MAJORS',TabNames{TabMelts}));
            % [A,ElementIndicies4Target] = ismember(targetStrings_Majors, PetrogenMajorElement_Strings);
            % nodata = find(ElementIndicies4Target==0);
            % ElementIndicies4Target(ElementIndicies4Target == 0) = max(ElementIndicies4Target);
            % ReorderedData = data2reorder(:,ElementIndicies4Target);
            % %pads rows of NaNs for elements with missing data
            % ReorderedData(:,nodata)=[NaN];
            %
            % MajorElements=ReorderedData(compIndicies,:);
            % InfoElement=InfoElement_ALL(compIndicies,infoIndicies);
            % TraceElements_ALL = eval(sprintf('%s_SM_TRACE',TabNames{TabMelts}));
            % TraceElements = TraceElements_ALL(compIndicies,:);
            %
            %     assignin('caller', [labelhere,'_SM'], sumMgNumNorm(MajorElements));
            %     assignin('caller', [labelhere,'_SM_Trace'],TraceElements);
            
            
%             assignin('caller',sprintf('Petrogen_TESource_%s',labelhere), Petrogen_Co);
%             assignin('caller',sprintf('Petrogen_TESourceName_%s',labelhere), Petrogen_Co_Name);
%             assignin('caller',sprintf('Petrogen_DataRatios_TESourceName_%s',labelhere), Petrogen_Co_Name);

            
            
            
            test= {labelhere [labelhere,'_Trace'] [labelhere,'_Info'] [labelhere,'_InfoALL'] [labelhere,'_Tormey'] [labelhere,'_DataRatios']...
                [labelhere,'_Residue_Trace'] [labelhere,'_SolidResidue_Trace'] [labelhere,'_CPX'] [labelhere,'_Ds']...
                [labelhere,'_Residue_DataRatios'] [labelhere,'_SolidResidue_DataRatios'] [labelhere,'_CPX_DataRatios'] ...
                ['DataRatios_', labelhere] [labelhere,'_DISEQ']...
                [sprintf('Petrogen_TESource_%s',labelhere)] [sprintf('Petrogen_TESourceName_%s',labelhere)] [sprintf('Petrogen_DataRatios_TESource_%s',labelhere)]... 
                };
            
            Names2Saves = {Names2Saves{:} test{:}};
            
            % [ExtraInfo] = reorderElements_Rows(SpreadsheetName, TabNames{i}, 'PlotOn', extrainfonames,'FirstInfo','LastInfo',column2unnormalize );
            % assignin('caller',sprintf('%s_Info',TabNames{i}),ExtraInfo);
            %
            
            
            %eval(sprintf('save(''%s'', ''%s'', ''%s'', ''%s'')', sprintf('%s',TabNames{i}), sprintf('%s',TabNames{i}),sprintf('%s_Tormey',TabNames{i}),sprintf('%s_Colors',TabNames{i})))
            All_INFO(kk,:) = InfoElement(1,:);
            
        end
    end
end

ElementRatios_Strings{21,1}='Eu';
ElementRatios_Strings{21,2}='Eu*';


% PetrogenTraceElement_Strings
% Petrogen_Co
% Petrogen_Co_Name
% Petrogen_BulkCompositionsNormalized
% Petrogen_Dmatrix
Save_PUM_WHDMM = {'TabName_Step2' 'ConstraintOptions' 'Info_Order' 'ElementRatiosIndex'  'ElementRatios_Strings' ...
      'PetrogenTraceElement_Strings' 'Petrogen_BulkCompositionsNormalized' 'Petrogen_Dmatrix' 'PetrogenMajorElement_Strings_Extracted'};

%Save_PUM_WHDMM = {'TabName_Step2' 'ConstraintOptions' 'Info_Order' 'ElementRatiosIndex','TraceElements4Normalization' 'DataRatios_TraceElements4Normalization' 'DataRatios_TraceDMM_TraceElements4Normalization' 'TraceElements4Normalization_TraceDMM' 'ElementRatios_Strings'};
Names2Saves = {Names2Saves{:} Save_PUM_WHDMM{:}};


save([pwd '/' subfolder '/' WorksheetName2Save '.mat'],Names2Saves{:})


fprintf('Name of extracted Step2 .mat file:<strong>''%s''</strong>',WorksheetName2Save)


end

